Developmental and Environmental Influences on Thalamocortical Connectivity are Patterned Along the S-A Axis in a Clinically Enriched Sample

Read in Data for Developmental and Environmental Characterization

Glasser parcel names/tracts

#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) 

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per bin

NeuroSynth cognitive atlas term maps

#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) 
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rh

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific tract list

tractlist.hbn <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HBN/HBN_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

Developmental statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir

#list all FA development files in development results directory for accessing
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F) 
files <- files[files %>% grep("hcpd", ., invert = T)] #only files for HBN and (to assess replication) PNC

filenames <- c()
#generate variable names to assign files data to 
for(name in files){
  Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
  filenames <- append(filenames, Rname) #save filenames into a character vector 
}

#read in files and assign to variables
for(i in 1:length(filenames)){
  
  Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
  
  if(grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.regions, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, S.A.axis, by = "tract")
    assign(Rfilename, x) 
  }
  
}
rm(x)

Environment statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/") #gam output dir, from gam_functions/fit_envGams.R

#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F) 
files <- files[files %>% grep("hbn", ., invert = F)] #only files for HBN and (to assess replication) PNC

for(i in 1:length(files)){
  filename <- files[i]
  Rname <- gsub('.{4}$', '', filename)
  x <- read.csv(filename, header = TRUE)
  assign(Rname, x) 
}
rm(x)

A Spectrum of Thalamocortical Developmental Trajectories

Cortex-wide thalamic connectivity developmental profiles

Number of significant developmental effects

#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

HBN

sig.effects(measure = "FA", atlas = "glasser", dataset = "hbn")
## [1] "In hbn, 166 thalamocortical connections (73.78 percent) show significant age-related changes in FA"

Thalamocortical connection GAM smooth functions

smoothcentered_FA_glasser_hbn.plot <- merge((smoothcentered_FA_glasser_hbn %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hbn %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_hbn.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.017, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(6, 8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Developmental_Trajectories_HBN.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.6)

Maturational patterns are highly reproducible

Correlation of thalamocortical connection development effect sizes between datasets

Pearson’s R

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hbn, by = c("tract", "label", "orig_parcelname"), suffixes = c("_pnc", "_hbn"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hbn)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hbn
## t = 15.616, df = 218, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6575492 0.7835828
## sample estimates:
##       cor 
## 0.7266226

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hbn, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 5e-05

Correlation plot

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hbn)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.051, 0.15)) +
  labs(x="\nPNC", y="HBN\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/FA_Ageeffect_correlation_PNCHBN.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.6)

Correlation of thalamocortical connection age of maturation between datasets

Pearson’s R

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hbn)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hbn
## t = 6.029, df = 151, p-value = 1.214e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3029697 0.5600095
## sample estimates:
##       cor 
## 0.4404722

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hbn, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0.0015

Early and Late Maturing Thalamocortical Connections

NeuroSynth decoding of age of maturation maps

Compute cognitive term-development map spatial correlations

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

Developmental timing decoding by Neurosynth cognitive terms

HBN

gameffects_neurosynth_hbn <- merge(gameffects_FA_glasser_hbn, neurosynth.terms, by = "label")

devpattern.neurosynth.hbn <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hbn, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))
devpattern.neurosynth.hbn <- rbind(slice_max(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.hbn
##                  term term.correlation
## 1   cognitive_control        0.6111711
## 2              recall        0.5162633
## 3             thought        0.5000909
## 4            decision        0.4826017
## 5          expectancy        0.4703547
## 6              memory        0.4511850
## 7           retrieval        0.4442939
## 8    memory_retrieval        0.4397963
## 9             context        0.4276347
## 10           judgment        0.4230770
## 11       multisensory       -0.5788485
## 12  visual_perception       -0.5575273
## 13 object_recognition       -0.5412269
## 14               gaze       -0.5204399
## 15       localization       -0.4933611
## 16          detection       -0.4488658
## 17         perception       -0.4268284
## 18  facial_expression       -0.3828248
## 19           movement       -0.3641246
## 20   visual_attention       -0.3538157
ggplot(devpattern.neurosynth.hbn, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 1.8, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.6,0.62), breaks = c(-0.6, -0.3, 0, 0.3, 0.6)) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_blank(),
   panel.grid.minor.y = element_blank(),
   panel.grid.minor.x = element_blank(),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Neurosynth_Maturationdecoding_HBN.pdf", device = "pdf", dpi = 500, width = 1.9, height = 3.1)

Overlapping terms across datasets

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))

devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map
merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hbn, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)
##                  term
## 1          expectancy
## 2   cognitive_control
## 3    memory_retrieval
## 4             thought
## 5              memory
## 6              recall
## 7   facial_expression
## 8          perception
## 9                gaze
## 10  visual_perception
## 11 object_recognition

Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis

gameffects_FA_glasser_hbn <- merge(gameffects_FA_glasser_hbn, (S.A.axis %>% select(tract, SA.axis)), by = "tract")

Thalamocortical connections mature at progressively older ages across the S-A axis

HBN

derivatives_FA_glasser_hbn$significant.derivative[derivatives_FA_glasser_hbn$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hbn, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(linewidth = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

Spearman’s rho

cor.test(gameffects_FA_glasser_hbn$SA.axis, gameffects_FA_glasser_hbn$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hbn$SA.axis and gameffects_FA_glasser_hbn$smooth.increase.offset
## S = 196170, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6899525

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.00205

Correlation plot

ggplot(gameffects_FA_glasser_hbn, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(11, 20), breaks = c(12, 14, 16, 18, 20)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/SAaxis_agematuration_hbn.pdf", device = "pdf", dpi = 500, width = 1.9, height = 1.6)

Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity

Plasticity marker development data

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)
#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)
EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 

EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)
df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts

HBN

plasticity.dev.hbn <- left_join(plasticity.dev, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$EI.ageslope
## S = 272148, p-value = 8.219e-15
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5698692

Spatial permutation (spin) based p-value

EI.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hbn, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$myelin.annualroc
## S = 998646, p-value = 2.643e-15
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.5783632

Spatial permutation (spin) based p-value

myelin.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hbn, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$amplitude.agedecline
## S = 147055, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6852288

Spatial permutation (spin) based p-value

amplitude.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

Correlation plot

ggplot(plasticity.dev.hbn, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.hbn, myelin.thalamus.hbn, amplitude.thalamus.hbn)
plasticity.ps
## [1] 0.00280 0.01495 0.00450
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.00675 0.01495 0.00675

Lolly plot of correlation strength

plasticity.results <- data.frame(map = factor(), corr = double())
plasticity.results <- plasticity.results %>% add_row(map = "E/I Ratio", corr = 0.5698692)
plasticity.results <- plasticity.results %>% add_row(map = "T1/T2 Ratio", corr = -0.5783632)
plasticity.results <- plasticity.results %>% add_row(map = "BOLD Amplitude", corr = 0.6852288)
plasticity.results$map <- factor(plasticity.results$map, ordered = TRUE, levels = c( "E/I Ratio", "T1/T2 Ratio", "BOLD Amplitude"))

ggplot(plasticity.results, aes(x = map, y = corr)) + 
  geom_segment(aes(x = map, xend = map,  yend =  
                     corr, y = 0), linewidth = 0.2, color = "#989898") +
  geom_point(size = 2, alpha = 0.85, color = "#FCAB6A", fill = "#FCAB6A", shape = 22) +
  labs(x="Cortical Development Map") +
  labs(y="\nThalamocortical Correlation (r)") +
  theme_classic() +
  scale_y_continuous(breaks = c(-0.6, -0.3, 0, 0.3, 0.6), limits = c(-0.62, 0.7)) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.text.x = element_text(vjust = 0.7, angle = 20),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Plasticitymap_correlation_HBN.pdf", device = "pdf", dpi = 500, width = 1.9, height = 1.58)

Associations between Neighborhood Socioeconomic Conditions and Thalamocortical Connectivity

Number of significant environmental effects

Significant neighborhood-level environment effects (neighborhood SES)

generalSES_maineffects_FA_glasser_hbn$significant <- p.adjust(generalSES_maineffects_FA_glasser_hbn$GAM.cov.pvalue, method = c("fdr")) < 0.05
sum(generalSES_maineffects_FA_glasser_hbn$significant)
## [1] 120
sum(generalSES_maineffects_FA_glasser_hbn$significant)/length(generalSES_maineffects_FA_glasser_hbn$significant)
## [1] 0.5333333

Extract significant connections

sigeffects.df <- generalSES_maineffects_FA_glasser_hbn %>% filter(significant == TRUE)

Percent of positive significant neighborhood-level environment effects

(sum(sigeffects.df$significant > 0)/nrow(sigeffects.df))*100
## [1] 100

Neighborhood Environment Effects Vary across the Sensorimotor-Association Axis

S-A axis correlation

Spearman’s rho

spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)

cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  spin.data$GAM.cov.tvalue and spin.data$SA.axis
## S = 199426, p-value = 0.0006705
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3075005

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.0349

Correlation plot

ggplot() +
  geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
  geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  labs(x="\nS-A axis", y="Environment effect (t)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/EnvEffect_SArank_correlationplot_HBN.pdf", device = "pdf", width = 1.9, height = 1.6)

Effect magnitude enrichment testing (S-A axis quintiles)

Mean environment effect t-value by S-A axis bin (empirical)

#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>% 
                             group_by(SA.axis.bin) %>% 
                             do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
                             unnest(cols = c("mean.tvalue.significant"))
envSES_effects_SAaxisbins
## # A tibble: 5 × 2
##   SA.axis.bin mean.tvalue.significant
##         <dbl>                   <dbl>
## 1           1                    2.84
## 2           2                    2.88
## 3           3                    2.92
## 4           4                    2.97
## 5           5                    3.22

Null distribution and enrichment testing by S-A axis quintile

#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)

## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix 
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full

## spin the empirical data 
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
  for (r in 1:nperm) {
    for (i in 1:nroi) {
      x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
      y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
    }
  }

## compute mean t-value for each S-A axis bin using spatially permuted x inputs 
### set up spun (x) and empirical (y) df
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm.empiricaly <- cbind(y, x.spatialperm) %>% as.data.frame()
colnames(x.spatialperm.empiricaly)[1] <- c("empirical.y")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedx <- sapply(x.spatialperm.empiricaly[2:ncol(x.spatialperm.empiricaly)], #for all permuted bins
                                              function(p){
                                                x.spatialperm.empiricaly %>% group_by(as.character(p)) %>% #group by permuted S-A axis bins
                                                do(mean.tvalue.significant = mean(.$empirical.y, na.rm = TRUE)) %>% #calculate mean t-value
                                                unnest(cols = c("mean.tvalue.significant")) %>% t()}) %>% 
                                                t() %>% as.data.frame()
                                                 
envSES_effects_SAaxisbins_permutedx <- envSES_effects_SAaxisbins_permutedx %>% mutate_if(is.character, as.numeric) #make it numeric
envSES_effects_SAaxisbins_permutedx <-  envSES_effects_SAaxisbins_permutedx[ , !c(TRUE,FALSE) ] 
colnames(envSES_effects_SAaxisbins_permutedx) <- c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue") 

## compute mean t-value for each S-A axis bin using spatially permuted y inputs 
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")

### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
                                       dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
                                       select(-empirical.x) %>% t() %>% as.data.frame() %>%
                                       set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))

histogram.data <- rbind(envSES_effects_SAaxisbins_permutedx, envSES_effects_SAaxisbins_permutedy)
saveRDS(histogram.data, "/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")
histogram.data <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")
envSES_effects_SAaxisbins_permutedx <- histogram.data[1:10000,]
envSES_effects_SAaxisbins_permutedy <- histogram.data[10001:20000,]
nperm = dim(perm.id.full)[2]
#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(quintile, enrichment_type){
  
  SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == quintile) %>% select(mean.tvalue.significant)
  
  #Histogram
  SA.colorvector <- c("#FFC228", "#F9D58C", "#e9dcf2", "#BE93C4", "#6F1382")
  
  null.tvalue.histogram <- histogram.data %>% select(sprintf("SA.bin%s.tvalue", quintile)) %>%
    ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[quintile]) +
    theme_classic() +
    geom_vline(xintercept =  SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
    xlab("Environment effect (t-value)") +
    ylab("Number of nulls") +
    scale_x_continuous(breaks = c(2, 3, 4), limits = c(2, 4)) +
    theme(
     legend.position = "none", 
     axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
     axis.text.y = element_blank(),
     axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
     axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
     axis.line.x = element_line(linewidth = .2), 
     axis.line.y = element_blank(),
     axis.ticks = element_line(linewidth = .2))
    
  #Permutation-based p
    if(enrichment_type == "smaller"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
    }
  
  if(enrichment_type == "larger"){
      p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- p.perm.xy
      p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
      perm.p <- (p.perm.xy+p.perm.yx)/2
  }
  
  permutation.results <- list(null.tvalue.histogram, perm.p)
  return(permutation.results)
}

S-A quintile 1

enveffects_SAbin_enrichment(quintile = 1, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.19655
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile1_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

S-A quintile 2

enveffects_SAbin_enrichment(quintile = 2, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.1787
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile2_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

S-A quintile 3

enveffects_SAbin_enrichment(quintile = 3, enrichment_type = "smaller")
## [[1]]

## 
## [[2]]
## [1] 0.19985
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile3_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

S-A quintile 4

enveffects_SAbin_enrichment(quintile = 4, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.5911
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile4_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)

S-A quintile 5

enveffects_SAbin_enrichment(quintile = 5, enrichment_type = "larger")
## [[1]]

## 
## [[2]]
## [1] 0.03145
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile5_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)